###Do Partisan Types Stop at the Water's Edge? 
#Partisan Types 1.R
#Josh Kertzer
#Last modified: February 13, 2020

### Note: this file runs the analysis in the main text - run Partisan Types 0.R first
### All of the following analyses were carried out using R version 3.5.2 on a 3 GHz Intel Xeon W iMac Pro running macOS Mojave 10.14.6, and replicated using R version 3.5.2 on a 2 GHz Intel Core i7 Macbook Pro running macOS Mojave 10.14.5.

###### Functions to generate stereotype measures #### 

#Stereotype content
s.content2 <- function(x, itemLevel=TRUE){
	temp <- x
	#First, re-order scale so that both, neither are the neutral midpoint
	temp[which(temp >=5)] <- 2.5 #Create new neutral midpoint
	temp[which(temp>=3)] <- temp[which(temp>=3)]+1 #Increase Republican responses by one
	temp[which(temp==2.5)] <- 3 #Now, move neutral midpoint up
	if (itemLevel==TRUE){
		return(temp)
	}
	else if(itemLevel==FALSE){
		return(mean(temp, na.rm=TRUE))
	}
}

#Stereotype intensity
s.intensity2 <- function(x, itemLevel=TRUE){
		temp <- x
		#First, re-order scale so that both, neither are the neutral midpoint
		temp[which(temp >=5)] <- 2.5 #Create new neutral midpoint
		temp[which(temp>=3)] <- temp[which(temp>=3)]+1 #Incerase Republican responses by one
		temp[which(temp==2.5)] <- 3 #Now, move neutral midpoint up
		#Now, calculate intensity
		if(itemLevel==TRUE){
			return(abs(temp - 3)/2)
		}
		else if(itemLevel == FALSE){
			return(mean(abs(temp-3)/2, na.rm=TRUE))
		}
}


#### Generate policy preference dataframes ####

##Study 1
#Note: Each column is for a different treatment
policy1 <- as.data.frame(matrix(NA, nrow=nrow(dat1), ncol=28))
for (i in 1:3){
	policy1[which(dat1$stopPiracy==(i-1)),i] <- dat1$Q4[which(dat1$stopPiracy==(i-1))] #Q4
	policy1[which(dat1$enviroSanction==(i-1)),i+3] <- dat1$Q5[which(dat1$enviroSanction==(i-1))]#Q5
	policy1[which(dat1$troopsOil==(i-1)),i+6] <- dat1$Q6[which(dat1$troopsOil==(i-1))]#Q6
	policy1[which(dat1$troopsCiv==(i-1)),i+9] <- dat1$Q7[which(dat1$troopsCiv==(i-1))]#Q7	
}
for (i in 1:2){
	policy1[which(dat1$freeTrade==(i-1)),i+12] <- dat1$Q8[which(dat1$freeTrade==(i-1))] #Q8
	policy1[which(dat1$proEnviro==(i-1)),i+14] <- dat1$Q9[which(dat1$proEnviro==(i-1))]#Q9
	policy1[which(dat1$interventionist==(i-1)),i+16] <- dat1$Q10[which(dat1$interventionist==(i-1))]#Q10
	policy1[which(dat1$armsControl==(i-1)),i+18] <- dat1$Q11[which(dat1$armsControl==(i-1))]#Q11	
	policy1[which(dat1$milSpend==(i-1)),i+20] <- dat1$Q12[which(dat1$milSpend==(i-1))]#Q12
	policy1[which(dat1$gunControl==(i-1)),i+22] <- dat1$Q13[which(dat1$gunControl==(i-1))]#Q13
	policy1[which(dat1$taxHike==(i-1)),i+24] <- dat1$Q14[which(dat1$taxHike==(i-1))]#Q14	
	policy1[which(dat1$abortions==(i-1)),i+26] <- dat1$Q15[which(dat1$abortions==(i-1))]#Q15	
}

policy1$id <- seq_len(nrow(policy1))

##Study 2
#Note: Each column is for a different treatment
policy2 <- as.data.frame(matrix(NA, nrow=nrow(dat2), ncol=22))
for (i in 1:3){
	policy2[which(dat2$troopsOil==(i-1)),i] <- dat2$Q157[which(dat2$troopsOil==(i-1))]#Q6
}
for (i in 1:2){
	policy2[which(dat2$freeTrade==(i-1)),i+3] <- dat2$Q161[which(dat2$freeTrade==(i-1))] 
	policy2[which(dat2$proEnviro==(i-1)),i+5] <- dat2$Q163[which(dat2$proEnviro==(i-1))]
	policy2[which(dat2$interventionist==(i-1)),i+7] <- dat2$Q165[which(dat2$interventionist==(i-1))]
	policy2[which(dat2$armsControl==(i-1)),i+9] <- dat2$Q167[which(dat2$armsControl==(i-1))]	
	policy2[which(dat2$milSpend==(i-1)),i+11] <- dat2$Q169[which(dat2$milSpend==(i-1))]
	policy2[which(dat2$gunControl==(i-1)),i+13] <- dat2$Q171[which(dat2$gunControl==(i-1))]
	policy2[which(dat2$taxHike==(i-1)),i+15] <- dat2$Q173[which(dat2$taxHike==(i-1))]	
	policy2[which(dat2$abortions==(i-1)),i+17] <- dat2$Q175[which(dat2$abortions==(i-1))]
	policy2[which(dat2$immig==(i-1)),i+19] <- dat2$Q180[which(dat2$immig==(i-1))]
	policy2[which(dat2$affAction==(i-1)),i+21] <- dat2$Q183[which(dat2$affAction==(i-1))]	
}

policy2$id <- seq_len(nrow(policy2))


#### Generate stereotype dataframes ####

## Study 1
#Note: Each column is for a different treatment
stereo1 <- as.data.frame(matrix(NA, nrow=nrow(dat1), ncol=28))
for (i in 1:3){
	stereo1[which(dat1$stopPiracy==(i-1)),i] <- dat1$Q17[which(dat1$stopPiracy==(i-1))] #Q17
	stereo1[which(dat1$enviroSanction==(i-1)),i+3] <- dat1$Q18[which(dat1$enviroSanction==(i-1))]#Q18
	stereo1[which(dat1$troopsOil==(i-1)),i+6] <- dat1$Q19[which(dat1$troopsOil==(i-1))]#Q19
	stereo1[which(dat1$troopsCiv==(i-1)),i+9] <- dat1$Q20[which(dat1$troopsCiv==(i-1))]#Q20
}
for (i in 1:2){
	stereo1[which(dat1$freeTrade==(i-1)),i+12] <- dat1$Q21[which(dat1$freeTrade==(i-1))] #Q21
	stereo1[which(dat1$proEnviro==(i-1)),i+14] <- dat1$Q22[which(dat1$proEnviro==(i-1))]#Q22
	stereo1[which(dat1$interventionist==(i-1)),i+16] <- dat1$Q23[which(dat1$interventionist==(i-1))]#Q23
	stereo1[which(dat1$armsControl==(i-1)),i+18] <- dat1$Q24[which(dat1$armsControl==(i-1))]#Q24
	stereo1[which(dat1$milSpend==(i-1)),i+20] <- dat1$Q25[which(dat1$milSpend==(i-1))]#Q25
	stereo1[which(dat1$gunControl==(i-1)),i+22] <- dat1$Q26[which(dat1$gunControl==(i-1))]#Q26
	stereo1[which(dat1$taxHike==(i-1)),i+24] <- dat1$Q27[which(dat1$taxHike==(i-1))]#Q27
	stereo1[which(dat1$abortions==(i-1)),i+26] <- dat1$Q28[which(dat1$abortions==(i-1))]#Q28	
}

## Study 2
#Note: each column is for a different treatment
stereo2 <- as.data.frame(matrix(NA, nrow=nrow(dat2), ncol=22))
for (i in 1:3){
	stereo2[which(dat2$troopsOil==(i-1)),i] <- dat2$Q183.1[which(dat2$troopsOil==(i-1))]
}
for (i in 1:2){
	stereo2[which(dat2$freeTrade==(i-1)),i+3] <- dat2$Q187[which(dat2$freeTrade==(i-1))]
	stereo2[which(dat2$proEnviro==(i-1)),i+5] <- dat2$Q189[which(dat2$proEnviro==(i-1))]
	stereo2[which(dat2$interventionist==(i-1)),i+7] <- dat2$Q191[which(dat2$interventionist==(i-1))]
	stereo2[which(dat2$armsControl==(i-1)),i+9] <- dat2$Q193[which(dat2$armsControl==(i-1))]
	stereo2[which(dat2$milSpend==(i-1)),i+11] <- dat2$Q195[which(dat2$milSpend==(i-1))]
	stereo2[which(dat2$gunControl==(i-1)),i+13] <- dat2$Q197[which(dat2$gunControl==(i-1))]
	stereo2[which(dat2$taxHike==(i-1)),i+15] <- dat2$Q199[which(dat2$taxHike==(i-1))]
	stereo2[which(dat2$abortions==(i-1)),i+17] <- dat2$Q201[which(dat2$abortions==(i-1))]
	stereo2[which(dat2$immig==(i-1)),i+19] <- dat2$Q181[which(dat2$immig==(i-1))]
	stereo2[which(dat2$affAction==(i-1)),i+21] <- dat2$Q182[which(dat2$affAction==(i-1))]	
}

##### Figure 1: Exploring different measures of partisan type
content1 <- rescale(as.data.frame(apply(stereo1,2,s.content2, itemLevel=TRUE)))
intensity1 <- rescale(as.data.frame(apply(stereo1,2,s.intensity2, itemLevel=TRUE)))

vector.m <- na.omit(melt(content1))
levels(vector.m$variable) <- c("Stop piracy", "Stop piracy (uni)", "Stop piracy (multi)", "Enviro sanction", "Enviro sanction (uni)", "Enviro sanction (multi)", "Troops oil", "Troops oil (uni)", "Troops oil (multi)", "Troops humanitarian", "Troops humanitarian (uni)", "Troops humanitarian (multi)", "Protectionist", "Free trade", "Anti-environment", "Pro-environment", "Isolationist", "Interventionist", "Anti-arms control", "Pro-arms control", "Decrease military spending", "Increase military spending", "Anti-gun control", "Pro-gun control", "Anti-tax hike", "Pro-tax hike", "Anti-abortion", "Pro-abortion")

demo.plot <- subset(vector.m, vector.m$variable=="Troops humanitarian" | vector.m$variable=="Anti-abortion" | vector.m$variable=="Pro-abortion")
levels(demo.plot$variable)[c(7,17,8)] <- c("Policy A", "Policy B", "Policy C")
demo.plot$variable <- relevel(demo.plot$variable, ref="Policy C")
demo.plot$variable <- relevel(demo.plot$variable, ref="Policy B")
demo.plot$variable <- relevel(demo.plot$variable, ref="Policy A")


content.vals <- content1[,c(27,10,28)]
intensity.vals <- intensity1[,c(27,10,28)]

B <- 1500
set.seed(43215)
plot.vals <- array(NA, dim=list(3,B,2)) #Variable, bootstraps, content/intensity
for (i in seq_len(B)){
	for (j in 1:3){ #Variables
		resample <- sample(1:nrow(content.vals), nrow(content.vals), replace=T)
		temp <- content.vals[resample,j]
		plot.vals[j,i,1] <- mean(temp,na.rm=TRUE)
		temp <- intensity.vals[resample,j]
		plot.vals[j,i,2] <- mean(temp,na.rm=TRUE)
	}
}

plot.vals2 <- data.frame(variable=c("Policy A", "Policy B", "Policy C"), mean=round(c(apply(plot.vals[,,1],1,mean,na.rm=TRUE), apply(plot.vals[,,2],1,mean,na.rm=TRUE)), digits=2), statistic=rep(c("Content", "Intensity"), each=3))


dev.new(height=3, width=7)
ggplot(demo.plot, aes(x=value)) + geom_bar(fill="darkgrey") + facet_wrap(~variable) + theme_bw() + labs(x="Distribution of responses by policy proposal", y="") + theme(axis.text.y=element_blank(), axis.ticks.y=element_blank()) + scale_x_continuous(breaks=c(0,0.5,1),labels=c("0"="0\nDefinitely\nDem", "0.5"="0.5\nBoth/\nNeither","1"="1\nDefinitely\nRep")) 

#Statistics in text and figure caption:
plot.vals2

##### Figures 2-3: Stereotype Content

## Figure 2: stereotype content in study 1
content1 <- rescale(as.data.frame(apply(stereo1,2,s.content2, itemLevel=TRUE)))
content1$id <- seq_len(nrow(content1))

B <- 1500
set.seed(43215)
content.boot1 <- array(NA, dim=list(ncol(content1)-1,B)) #Variable, bootstraps
for (i in seq_len(B)){
	for (j in 1:28){ #Variables
		resample <- sample(content1$id, nrow(content1), replace=T)
		temp <- content1[resample,j]
		content.boot1[j,i] <- mean(temp,na.rm=TRUE)
	}
}

content.mat1 <- data.frame(V1=apply(content.boot1,1,mean,na.rm=TRUE), V2=apply(content.boot1,1,quantile,0.025,na.rm=TRUE), V3=apply(content.boot1,1,quantile,0.975,na.rm=TRUE), V4=c(1:28))
content.mat1$V5 <- "Domestic"
content.mat1$V5[c(1:12,13:14,17:20,21:22)] <- "Foreign" 
content.mat1$V6 <- factor(c("Stop piracy", "Stop piracy (unilateral)", "Stop piracy (multilateral)", "Enviro sanction", "Enviro sanction (unilateral)", "Enviro sanction (multilateral)", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)", "Troops humanitarian", "Troops humanitarian (unilateral)", "Troops humanitarian (multilateral)", "Protectionist", "Free trade", "Anti-environment", "Pro-environment", "Isolationist", "Interventionist", "Anti-arms control", "Pro-arms control", "Decrease military spending", "Increase military spending", "Anti-gun control", "Pro-gun control", "Anti-tax hike", "Pro-tax hike", "Anti-abortion", "Pro-abortion"), levels=c("Pro-abortion", "Anti-abortion", "Pro-tax hike", "Anti-tax hike", "Pro-environment", "Anti-environment", "Pro-gun control", "Anti-gun control", "Pro-arms control", "Anti-arms control", "Interventionist", "Isolationist", "Free trade", "Protectionist", "Increase military spending", "Decrease military spending", "Troops humanitarian", "Troops humanitarian (unilateral)", "Troops humanitarian (multilateral)", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)", "Enviro sanction", "Enviro sanction (unilateral)", "Enviro sanction (multilateral)", "Stop piracy", "Stop piracy (unilateral)", "Stop piracy (multilateral)"))
content.mat1$V7 <- content.mat1$V5
content.mat1$V7[1:12] <- "Foreign (Approach)"
content.mat1$V8 <- c(rep(1:4,each=3), rep(5:12,each=2))

dev.new(height=7, width=8.8)
ggplot(data=content.mat1, aes(x=V1,y=V6)) + geom_point() + geom_segment(aes(x=V2, y=V6,xend=V3, yend=V6)) +theme_bw() + guides(color=FALSE, alpha=FALSE) + scale_color_grey() + facet_wrap(~V7, nrow=3, scales="free") +  labs(x="Stereotype content", y="") + geom_line(aes(x=V1, y=V6, group=V8, alpha=0.8)) + xlim(0,1)
#In post-production, edit the ordering of the connective lines in the bottom panel manually; add Democratic / Neutral / Republican legend on x axis

## Figure 3: stereotype content in study 2

content2 <- rescale(as.data.frame(apply(stereo2,2,s.content2, itemLevel=TRUE)))
content2$id <- seq_len(nrow(content2))

B <- 1500
set.seed(43215)
content.boot2 <- array(NA, dim=list(ncol(content2)-1,B)) #Variable, bootstraps
for (i in seq_len(B)){
	for (j in 1:23){ #Variables
		resample <- sample(content2$id, nrow(content2), replace=T)
		temp <- content2[resample,j]
		content.boot2[j,i] <- mean(temp,na.rm=TRUE)
	}
}

content.mat2 <- data.frame(V1=apply(content.boot2,1,mean,na.rm=TRUE), V2=apply(content.boot2,1,quantile,0.025,na.rm=TRUE), V3=apply(content.boot2,1,quantile,0.975,na.rm=TRUE), V4=c(1:23))
content.mat2$V5 <- "Domestic"
content.mat2$V5[c(1:3,4:5,8:13,22:23)] <- "Foreign" 
content.mat2$V6 <- factor(c("Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)","Protectionist", "Free trade", "Anti-environment", "Pro-environment", "Isolationist", "Interventionist", "Anti-arms control", "Pro-arms control", "Decrease military spending", "Increase military spending", "Anti-gun control", "Pro-gun control", "Anti-tax hike", "Pro-tax hike", "Anti-abortion", "Pro-abortion", "Anti-affirmative action", "Pro-affirmative action", "Anti-immigration", "Pro-immigration"), levels=c("Pro-affirmative action", "Anti-affirmative action", "Pro-abortion", "Anti-abortion", "Pro-tax hike", "Anti-tax hike", "Pro-environment", "Anti-environment", "Pro-gun control", "Anti-gun control", "Pro-immigration", "Anti-immigration", "Pro-arms control", "Anti-arms control", "Interventionist", "Isolationist", "Free trade", "Protectionist", "Increase military spending", "Decrease military spending", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)"))
content.mat2$V7 <- content.mat2$V5
content.mat2$V7[1:3] <- "Foreign (Approach)"
content.mat2$V8 <- c(rep(1,3), rep(2:11,each=2))

dev.new(height=7, width=8.8)
ggplot(data=content.mat2, aes(x=V1,y=V6)) + geom_point() + geom_segment(aes(x=V2, y=V6,xend=V3, yend=V6)) +theme_bw() + guides(color=FALSE, alpha=FALSE) + scale_color_grey() + facet_wrap(~V7, nrow=3, scales="free") +  labs(x="Stereotype content", y="") + geom_line(aes(x=V1, y=V6, group=V8, alpha=0.8)) + xlim(0,1)
#In post-production, edit the ordering of the connective lines in the bottom panel manually; add Democratic / Neutral / Republican legend on x axis

### Interpretation of results in text

#Abortion
round(content.mat1$V1[28], digits=2) #Pro-choice, study 1: 0.25
round(content.mat1$V1[27], digits=2) #Anti-abortion, study 1: 0.71
round(content.mat1$V1[28] - content.mat1$V1[27], digits=2) #-46 percentage point change in 2014
round(content.mat2$V1[19], digits=2) #Pro-choice, study 2: 0.27
round(content.mat2$V1[18], digits=2) #Pro-choice, study 2: 0.72
round(content.mat2$V1[19] - content.mat2$V1[18], digits=2) #Study 2: -44 percentage point change in 2018

#Average size of treatment effects for domestic issues, 2014: 40%
round((abs(content.mat1$V1[16] - content.mat1$V1[15]) + abs(content.mat1$V1[24] - content.mat1$V1[23]) + abs(content.mat1$V1[26] - content.mat1$V1[25]) + abs(content.mat1$V1[28] - content.mat1$V1[27]))/4, digits=2) #Average: 40%
#Average size of treatment effects for domestic issues, 2018: 37%
round((abs(content.mat2$V1[7] - content.mat2$V1[6]) + abs(content.mat2$V1[15] - content.mat2$V1[14]) + abs(content.mat2$V1[17] - content.mat2$V1[16]) + abs(content.mat2$V1[19] - content.mat2$V1[18]) + abs(content.mat2$V1[21] - content.mat2$V1[20]))/5, digits=2) #Average: 37%

#Immigration
round(content.mat2$V1[23] - content.mat2$V1[22], digits=2) #Pro immigration vs anti: -32%

#Military spending
round(content.mat1$V1[22] - content.mat1$V1[21], digits=2) #Increase military spending vs decrease in 2014: 30%
round(content.mat2$V1[13] - content.mat2$V1[12], digits=2) #Increase mil spending vs decrease in 2018: 24%

#Arms control
round(content.mat1$V1[20] - content.mat1$V1[19], digits=2) #Pro arms-control vs anti-arms control in 2014: -14%
round(content.mat2$V1[11] - content.mat2$V1[10], digits=2) #Pro arms-control vs anti-arms control in 2018: -9%

#Trade
round(content.mat1$V1[14] - content.mat1$V1[13], digits=2) #Free trade vs protectionism in 2014: -3%
round(content.mat2$V1[5] - content.mat2$V1[4], digits=2) #Free trade vs protectionism in 2018: -13%

#Isolationism
round(content.mat1$V1[18] - content.mat1$V1[17], digits=2) #Interventionist vs isolationist in 2014: 1%
round(content.mat2$V1[9] - content.mat2$V1[8], digits=2) #Interventionist vs isolationist in 2018: -10%

#Foreign policy approaches
round((abs(content.mat1$V1[3] - content.mat1$V1[2]) + abs(content.mat1$V1[6] - content.mat1$V1[5]) + abs(content.mat1$V1[9] - content.mat1$V1[8]) + abs(content.mat1$V1[12] - content.mat1$V1[11]))/4, digits=2) #Average effect in 2014: 3%
round(content.mat2$V1[3] - content.mat2$V1[2], digits=2) #Average effect in 2018: 3%

##### Figure 4: Stereotype intensity

## Figure 4a: Study 1
intensity1 <- rescale(as.data.frame(apply(stereo1,2,s.intensity2, itemLevel=TRUE)))
intensity1$id <- seq_len(nrow(intensity1))

B <- 1500
set.seed(43215)
intensity.boot1 <- array(NA, dim=list(ncol(intensity1)-1,B)) #Variable, bootstraps
for (i in seq_len(B)){
	for (j in 1:28){ #Variables
		resample <- sample(intensity1$id, nrow(intensity1), replace=T)
		temp <- intensity1[resample,j]
		intensity.boot1[j,i] <- mean(temp,na.rm=TRUE)
	}
}

intensity.mat1 <- data.frame(V1=apply(intensity.boot1,1,mean,na.rm=TRUE), V2=apply(intensity.boot1,1,quantile,0.025,na.rm=TRUE), V3=apply(intensity.boot1,1,quantile,0.975,na.rm=TRUE), V4=c(1:28))
intensity.mat1$V5 <- "Domestic"
intensity.mat1$V5[c(1:12,13:14,17:20,21:22)] <- "Foreign" 
intensity.mat1$V6 <- factor(c("Stop piracy", "Stop piracy (unilateral)", "Stop piracy (multilateral)", "Enviro sanction", "Enviro sanction (unilateral)", "Enviro sanction (multilateral)", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)", "Troops humanitarian", "Troops humanitarian (unilateral)", "Troops humanitarian (multilateral)", "Protectionist", "Free trade", "Anti-environment", "Pro-environment", "Isolationist", "Interventionist", "Anti-arms control", "Pro-arms control", "Decrease military spending", "Increase military spending", "Anti-gun control", "Pro-gun control", "Anti-tax hike", "Pro-tax hike", "Anti-abortion", "Pro-abortion"), levels=c("Pro-abortion", "Anti-abortion", "Pro-tax hike", "Anti-tax hike", "Pro-environment", "Anti-environment", "Pro-gun control", "Anti-gun control", "Pro-arms control", "Anti-arms control", "Interventionist", "Isolationist", "Free trade", "Protectionist", "Increase military spending", "Decrease military spending", "Troops humanitarian", "Troops humanitarian (unilateral)", "Troops humanitarian (multilateral)", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)", "Enviro sanction", "Enviro sanction (unilateral)", "Enviro sanction (multilateral)", "Stop piracy", "Stop piracy (unilateral)", "Stop piracy (multilateral)"))
intensity.mat1$V7 <- intensity.mat1$V5
intensity.mat1$V7[1:12] <- "Foreign (Approach)"
intensity.mat1$V8 <- c(rep(1:4,each=3), rep(5:12,each=2))

dev.new(height=5.2, width=8)
ggplot(data=intensity.mat1, aes(x=V1,y=V6)) + geom_point() + geom_segment(aes(x=V2, y=V6,xend=V3, yend=V6)) +theme_bw() + guides(color=FALSE, alpha=FALSE) + scale_color_grey() + facet_grid(V7 ~ 1, scales="free", space="free") +  labs(x="Stereotype intensity", y="") + xlim(0,1) #Change label of row in post-production


## Figure 4b: Study 2
intensity2 <- rescale(as.data.frame(apply(stereo2,2,s.intensity2, itemLevel=TRUE)))
intensity2$id <- seq_len(nrow(intensity2))

B <- 1500
set.seed(43215)
intensity.boot2 <- array(NA, dim=list(ncol(intensity2)-1,B)) #Variable, bootstraps
for (i in seq_len(B)){
	for (j in 1:23){ #Variables
		resample <- sample(intensity2$id, nrow(intensity2), replace=T)
		temp <- intensity2[resample,j]
		intensity.boot2[j,i] <- mean(temp,na.rm=TRUE)
	}
}

intensity.mat2 <- data.frame(V1=apply(intensity.boot2,1,mean,na.rm=TRUE), V2=apply(intensity.boot2,1,quantile,0.025,na.rm=TRUE), V3=apply(intensity.boot2,1,quantile,0.975,na.rm=TRUE), V4=c(1:23))
intensity.mat2$V5 <- "Domestic"
intensity.mat2$V5[c(1:3,4:5,8:13,22:23)] <- "Foreign" 
intensity.mat2$V6 <- factor(c("Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)","Protectionist", "Free trade", "Anti-environment", "Pro-environment", "Isolationist", "Interventionist", "Anti-arms control", "Pro-arms control", "Decrease military spending", "Increase military spending", "Anti-gun control", "Pro-gun control", "Anti-tax hike", "Pro-tax hike", "Anti-abortion", "Pro-abortion", "Anti-affirmative action", "Pro-affirmative action", "Anti-immigration", "Pro-immigration"), levels=c("Pro-affirmative action", "Anti-affirmative action", "Pro-abortion", "Anti-abortion", "Pro-tax hike", "Anti-tax hike", "Pro-environment", "Anti-environment", "Pro-gun control", "Anti-gun control", "Pro-immigration", "Anti-immigration", "Pro-arms control", "Anti-arms control", "Interventionist", "Isolationist", "Free trade", "Protectionist", "Increase military spending", "Decrease military spending", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)"))
intensity.mat2$V7 <- intensity.mat2$V5
intensity.mat2$V7[1:3] <- "Foreign (Approach)"
intensity.mat2$V8 <- c(rep(1,3), rep(2:11,each=2))

dev.new(height=5.2, width=7.62)
ggplot(data=intensity.mat2, aes(x=V1,y=V6)) + geom_point() + geom_segment(aes(x=V2, y=V6,xend=V3, yend=V6)) +theme_bw() + guides(color=FALSE, alpha=FALSE) + scale_color_grey() + facet_grid(V7 ~ 1, scales="free", space="free") +  labs(x="Stereotype intensity", y="") + xlim(0,1)
#Add intensity legend on x axis, combine two panels in post-production because of JOP's requirement that multi-panel figures be consolidated

##### Assessing the accuracy of partisan types

### Reshape data in order to estimate linear mixed models

#First, create mergeable dataset of policy preferences
policy1.m <- cbind(policy1[,1:28], V29=NA, V30=NA, V31=NA, V32=NA, id=policy1[,29]) #Add four new columns for the two extra policies (two stances for each)
policy2.m <- cbind(NA,NA,NA,NA,NA,NA,policy2[,1:3],NA,NA,NA, policy2[,4:23])
colnames(policy2.m)[1:32] <- paste("V",1:32,sep="")
policy2.m$id <- 1017:2021
policies <- rbind(policy1.m, policy2.m)
policies <- melt(policies, id.vars="id")

#Now, mergeable dataset of stereotype intensity
stereo1.m <- cbind(intensity1[,1:28], V29=NA, V30=NA, V31=NA, V32=NA, id=intensity1[,29]) #Add four new columns for the two extra policies (two stances for each)
stereo2.m <- cbind(NA,NA,NA,NA,NA,NA,intensity2[,1:3],NA,NA,NA, intensity2[,4:23])
colnames(stereo2.m)[1:32] <- paste("V",1:32,sep="")
stereo2.m$id <- 1017:2021

#Now, merge together
types <- rbind(stereo1.m, stereo2.m)
types <- melt(types, id.vars="id")
types <- cbind(types, policies$value)
colnames(types) <- c("id", "issue", "value", "prefs")

#Now add other covariates
types$year <- ifelse(types$id <= 1016,"2014", "2018")
types$prefsStr <- abs(types$prefs-4)/3 #Standardized measure of respondent's intensity of preferences on given issue

#Foreign policy issue?
foreign <- data.frame(t(rep(0,32)))
foreign[c(1:12,13:14,17:20,21:22,31:32)] <- 1 #Foreign policy dummy
colnames(foreign) <- paste("V",1:32,sep="")

#Polarization score: calculate difference between Democrats and Republicans on issue preferences
policy1.m.dems <- policy1.m[which(dat1$pid3==1),-33]
policy1.m.reps <- policy1.m[which(dat1$pid3==3),-33]
policy2.m.dems <- policy2.m[which(dat2$pid3==1),-33]
policy2.m.reps <- policy2.m[which(dat2$pid3==3),-33]

pol.1 <- abs(apply(policy1.m.reps,2,mean,na.rm=TRUE) - apply(policy1.m.dems,2,mean,na.rm=TRUE))/(7-1)
pol.2 <- abs(apply(policy2.m.reps,2,mean,na.rm=TRUE) - apply(policy2.m.dems,2,mean,na.rm=TRUE))/(7-1)

policies1 <- cbind(melt(foreign), melt(pol.1))
colnames(policies1) <- c("issue","Foreign", "Polarization")
policies2 <- cbind(melt(foreign), melt(pol.2))
colnames(policies2) <- c("issue","Foreign", "Polarization")
policies1$year <- "2014"
policies2$year <- "2018"
policies <- rbind(policies1, policies2)

#Merge foreign policy issue and polarization scores into the dataset
types <- merge(types, policies, by=c("issue","year"))

#Now, merge in other individual differences
loc.a <- match(c("male", "age", "educ1", "pid1", "ideo1", "white", "polInterest1"), colnames(dat1))
ind.diff1 <- dat1[,loc.a]
ind.diff1$year <- "2014"
ind.diff1$id <- seq_len(nrow(ind.diff1))
loc.b <- match(c("male", "age", "educ1", "pid1", "ideo1", "white", "polInterest1"), colnames(dat2))
ind.diff2 <- dat2[,loc.b]
ind.diff2$year <- "2018"
ind.diff2$id <- 1017:2021
ind.diff <- rbind(ind.diff1, ind.diff2)

types <- merge(types, ind.diff, by=c("id","year"))

### Table 1: Linear mixed models

mod.1 <- lmer(value ~ 1 + (1|id) + (1|issue) + (1|year), data=types) #More variance from issues than people
v <- VarCorr(mod.1)
s2 <- attr(v, "sc")^2
t0 <- as.numeric(v$id)
t1 <- as.numeric(v$issue)
t2 <- as.numeric(v$year)
t0/(s2+t0+t1+t2) #32.8% from respondent
t1/(s2+t0+t1+t2) #7.2% from issue
t2/(s2+t0+t1+t2) #0% from year
t0/t1 #4.55 times more variation across respondents than issues (reported in text)

mod.2 <- lmer(value ~ 1 + age + male + white + age + educ1 + pid1 + polInterest1 + (1|id) + (1|issue) + (1|year), data=types)

mod.3 <- lmer(value ~ 1 + age + male + white + age + educ1 + pid1 + polInterest1 + prefsStr + (1|id) + (1|issue) + (1|year), data=types)

mod.4 <- lmer(value ~ 1 + age + male + white + age + educ1 + pid1 + polInterest1 + Foreign + (1|id) + (1|issue) + (1|year), data=types) 

mod.5 <- lmer(value ~ 1 + age + male + white + age + educ1 + pid1 + polInterest1 + Polarization + (1|id) + (1|issue) + (1|year), data=types) 

stargazer(mod.1, mod.2, mod.3, mod.4, mod.5, title="Mixed linear model: stereotype strength", omit.stat=c("LL", "ser", "f"), style="apsr", digits=3, label="tab:stereotypicality", covariate.labels=c("Age", "Male", "White", "Education", "Partisanship", "Political interest", "Strength of preferences", "Foreign policy issue", "Polarization", "Constant"))

#What's the effect size of polarization?
summary(types$Polarization) #Min: 0, Max: 0.345
round(max(types$Polarization, na.rm=TRUE)*fixef(mod.5)[8], digits=2) #16%

#### Figure 5: Changes in second-order beliefs from 2014-18 reflect changes in first-order preferences
temp1 <- content1[,c(7:9,13:28)]
temp2 <- content2[,c(1:19)]
colnames(temp1) <- paste("V",1:19,sep="") #Need identical column names in order to rbind

stereo.dat <- as.data.frame(rbind(temp1, temp2))
stereo.dat$year <- c(rep(2014,nrow(content1)), rep(2018, nrow(content2)))
stereo.dat$id <- seq_len(nrow(stereo.dat))
stereo.dat$weight=c(rep(NA, nrow(content1)), dat2$weight)

temp1 <- policy1[,c(7:9,13:28)]
temp2 <- policy2[,c(1:19)]
colnames(temp1) <- paste("V",1:19,sep="") #Need identical column names in order to rbind

policy.dat <- as.data.frame(rbind(temp1, temp2))
policy.dat$year <- c(rep(2014,nrow(policy1)), rep(2018, nrow(policy2)))
policy.dat$id <- seq_len(nrow(policy.dat))
policy.dat$pid3 <- c(dat1$pid3, dat2$pid3)
policy.dat$weight=c(rep(NA, nrow(policy1)), dat2$weight)


dif.content <- rep(NA,19)
for (j in 1:19){ #Variables
	temp.18 <- which(stereo.dat$year==2018)
	temp.14 <- which(stereo.dat$year==2014)
	dif.content[j] <- mean(stereo.dat[temp.18,j],na.rm=TRUE) - mean(stereo.dat[temp.14,j],na.rm=TRUE)
	}

dif.policy <- rep(NA,19)
for (j in 1:19){ #Variables
	temp.18.d <- which(policy.dat$year==2018 & policy.dat$pid3==1)
	temp.14.d <- which(policy.dat$year==2014 & policy.dat$pid3==1)
	temp.18.r <- which(policy.dat$year==2018 & policy.dat$pid3==3)
	temp.14.r <- which(policy.dat$year==2014 & policy.dat$pid3==3)
	dif.policy[j] <- (mean(policy.dat[temp.18.r,j],na.rm=TRUE) - mean(policy.dat[temp.14.r,j],na.rm=TRUE)) - (mean(policy.dat[temp.18.d,j],na.rm=TRUE) - mean(policy.dat[temp.14.d,j],na.rm=TRUE))
}


mod.r <- lm(dif.content ~ dif.policy)
summary(mod.r)$r.squared #r2= 0.72	
cor(dif.content, dif.policy, use="complete.obs")

dif.plot <- data.frame(stereo=dif.content, policy=dif.policy, V6=factor(c("Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)", "Protectionist", "Free trade", "Anti-environment", "Pro-environment", "Isolationist", "Interventionist", "Anti-arms control", "Pro-arms control", "Decrease military spending", "Increase military spending", "Anti-gun control", "Pro-gun control", "Anti-tax hike", "Pro-tax hike", "Anti-abortion", "Pro-abortion"), levels=c("Pro-abortion", "Anti-abortion", "Pro-tax hike", "Anti-tax hike", "Pro-environment", "Anti-environment", "Pro-gun control", "Anti-gun control", "Pro-arms control", "Anti-arms control", "Interventionist", "Isolationist", "Free trade", "Protectionist", "Increase military spending", "Decrease military spending", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)"))) 

dev.new(height=7, width=7)
ggplot(dif.plot, aes(x=policy, y=stereo)) + geom_point() + theme_bw() + guides(alpha=FALSE) + labs(x="Difference-in-difference in policy preferences: Republicans vs Democrats, 2014-2018", y="Change in stereotype content, 2014-18") + geom_vline(xintercept=0, lty=3) + geom_hline(yintercept=0, lty=3) + geom_smooth() + geom_label_repel(aes(policy,stereo,label=V6), data=dif.plot, force=1, show.legend=FALSE, size=3) 
